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ABSTRACT 

We describe the CRASH (Center for Radiative Shock Hydrodynamics) code, a block adaptive 
mesh code for muhi-material radiation hydrodynamics. The implementation solves the radiation 
diffusion model with the gray or multigroup method and uses a flux limited diffusion approxi- 
mation to recover the free-streaming limit. The electrons and ions are allowed to have different 
temperatures and we include a flux limited electron heat conduction. The radiation hydro- 
dynamic equations are solved in the Eulerian frame by means of a conservative finite volume 
discretization in either one, two, or three-dimensional slab geometry or in two-dimensional cylin- 
drical symmetry. An operator split method is used to solve these equations in three substeps: 
(f ) solve the hydrodynamic equations with shock-capturing schemes, (2) a linear advection of the 
radiation in frequency-logarithm space, and (3) an implicit solve of the stiff radiation diffusion, 
heat conduction, and energy exchange. We present a suite of verification test problems to demon- 
strate the accuracy and performance of the algorithms. The CRASH code is an extension of the 
Block-Adaptive Tree Solarwind Roe Upwind Scheme (BATS-R-US) code with this new radiation 
transfer and heat conduction library and equation-of-state and multigroup opacity solvers. Both 
CRASH and BATS-R-US are part of the publicly available Space Weather Modeling Framework 
(SWMF). 

Subject headings: hydrodynamics - methods: numerical - radiative transfer 



INTRODUCTION 



As photons travel through matter, th e radiation field exp e rience s changes due to ne t total emission , 
absorption, and scattering, see for instance iMihalas fc Mihalas I (| 19841 ): iPomraming I (|2005f ) : iDrake I (|2006f ). 
At high enough energy density the radiation will heat and accelerate the plasma. This coupled system is 
radiation hydrodynamics. The radiation can at a fundamental level be described by the time evolution of 
the spectral radiation intensity Ii,{r, t, n, z^), which is the radiation energy per unit area, per unit solid angle 
in the direction of photon propagation n, per unit interval of photon frequency and per unit time interval. 
Several method have been developed to solve the radiation field in various degree of physics fidelity. 

In Monte Carlo radiative transfer methods, the radiation is statistically evaluated. Small photon pack- 
ets are created with their energy and propagation direction statistically selected. The packets are prop- 
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agated through n iatter using the radiation transfer equation (jNavakshin et al 
Baek et al. 1 120091 ). Characteristic methods use integration along rays of various lengths to solve for the 
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angular structure of the radiation transport. A recent conservative, causal ray-tracing method was devel- 
oped and combined wi th a short characteristic ray-tracing for the transfer calculations of ionizing radiation 



(jMellema et al. 1 120061 ) . Solar surface magneto-convection simulations are increasingly realistic and use a 



three-dimensional, non-gray, approximate local thermodynamic equilibrium (LTE), radiative transfer for the 
heating and cooling of plasm a. These simulations are typically formulated for four frequency b ins in the 



radiative transport equation (jVogler et al. 1 120051 : IStein et al. Il2007t iMartmez-Svkora et al 



20091) . 



For some applications, simplifications to the radiation transfer can be made by calculating moments 
of the radiation intensity over the solid angle 57. The spectral radiation energy and the spectral radiation 
energy flux are defined by the 0th and 1st moments as 



E^{r,t,iy) - j I^{r,t,n,iy)dil, F^{r, t,i^)= j nl^{r,t,n,i^)dn. 

In addition, the spectral radiation pressure tensor is defined by the second moment 

Py(r,t,i') = - / nnIi,[T,t,n,v)dVL. 

A whole class of radiation transfer models are based on solving the corresponding radiati on moment equa- 
tions , using a closure relation between the spectral pressure tensor and the spectral intensity (jMihalas fc Mihalas 
1984t IPomraming llioosi brake Ibooel ). 



(1) 



(2) 



A radiation-hydrodynamics code based on variable Eddington tensor (VET) methods (jStone et al. 



19921 ) can still capture the angular structure of the radiation field by relating the spectral radiation pressure 
tensor to the spectral radiation intensity and the method is applicable for both the optically thin and thick 
regime. Optically thin vers ions of the VET method have been used in the context of cosmological reionization 
|Petkova fc Springel lliooj ). 



Further simplification assumes that the radiation pressure is isotropic and proportional to the radiation 
energy. T his is the diffus ion approximation. Several codes have been developed using this approximation. 
HYDRA (jMarinak 1 1200 II ) is an arbitrary Lagrange Eulerian code for 2D and 3D radiation hydrodynamics. 
The radiation transfer model is based upon either flux limited multigroup or implicit Monte Carlo radiation 
transport. The Eulerian code RAGE (jGittings et al. 112008) uses a cell- based adaptive mesh refinement to 
achieve resolved radiative hydrodynamic flows. HYADE S (Larsen fc Lan e 1994) solve the radiation hydro- 
dynamic equations on a Lagrangian mesh, while GALE ( Barton 19851 ) can use either a fixed Eulerian mesh, 
an embedded Lagrangian mesh, or a partially embedded, partially remapped mesh. Our newly developed 
radiation-hydrodynamic solver uses an Eulerian grid together with a block-based adaptive mesh refinement 
strategy. 

We limit the discussion of the radiation hydrodynamics implementation in GRASH to plasmas in the 
absence of magnetic field. Most of the description in this paper can, however, easily be extended to magne- 
tohydrodynamic (MHD) plasmas as well. Indeed, since the GRA SH code is essentially the magnetohydrody- 



namic BATS-R-US code (jPowell et al. Ill999t iToth et al. 1120101 ) extended with libraries containing radiation 



transport, equation-of-state (EOS), and opacity solvers, the implementation for the coupling between the 
radiation field and MHD pla smas is read i ly ava ilable. The GRASH code uses the recently developed block 
adaptive tree library (BATL. lToth et al. I (j2010f )l. Here we will focus on the radiation implementation. Both 
the CRA SH and BATS-R-U S codes are publicly available as part of the Space Weather Modeling Framework 
(SWMF. iToth et al. I (|2005|)) or can be used as stand-alone codes. 



In the following. Section [2] introduces the radiation hydrodynamic equations for multi-material plasmas, 
in a form general enough to apply at high energy density. Section [3] describes the numerical algorithms 
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to solve these equations. Next, in Section |4] verification tests, are presented for radiation and electron heat 
conduction on non-uniform meshes in ID, 2D, and 3D slab geometry and in axially symmetric (rz) geometry. 
We also show a full system multi-material radiation hydrodynamic simulation on an adaptively refined mesh 
and demonstrate good scaling up to 1000 processors. The paper is summarized in Section O 



2. EQUATIONS OF RADIATION HYDRODYNAMICS IN DENSE PLASMAS 

The equations of radiation hydrodynamics describe the time evolution of both matter and radiation. 
For the applications that supported the work reported here, the code must be able to model matter as a high 
energy density plasma that is in LT E so that the population of all atomic and ion states can be obtained from 



statistical physics (see for instance iLandau fc Lifshitz I (jlOSOj)). We allow for multiple materials throughout 



the spatial domain of interest, but restrict the analysis to plasma flows that are far from relativistic. The 
materials can be heated to sufficiently high temperatures so that they can ionize and create free electrons, 
introducing the need for a time evolution equation for the electron energy density. The electrons transfer 
heat by thermal heat conduction and emit and absorb photon radiation. The radiation model discussed in 
this paper is non-equilibrium diffusion, in which the electron and radiation temperature can be different. 
We approximate the radiation transfer with a gray or multigroup fiux limited diffusion (FLD). This model 
is also of interest for application to a number of astrophysical problems. 

In the following subsections, we will describe the radiative transfer equations for the evolution of the 
multigroup radiation energy densities f Section l2.ip in the FLD approximation (Section l2.5p . The coupling of 
the radiation field to the two species hydrodynamic equations of electrons and ions are discussed in Section 
12.21 In Section [2731 the method for tracking the different materials is treated, while the lookup tables used 
for of the EOS and opacities are mentioned in Section [2741 



2.1. Radiation Transport 



In this section, we will build up the form of the radiation transport in the multigroup diffusion approxi- 
mation that is used for the implem entation in the CRASH co de. The spectral pressure tensor, equation 
is often approximated in the form (jMihalas fc Mihalas II1984I) 



where 



T,(r,i,iy) - i(l-x.)I 



1 F F 

2(3X.-1)^, 



(3) 



(4) 



is the spectral Eddington tensor, Xf is the Eddington factor, and I is the identity matrix. The second term 
on the right hand side is a dyad constructed from the direction of the spectral radiation flux. The pressure 
tensor can be used to arrive at a time evolution equation for the solid angle integrated spectral radiation 
energy (jBuchler lllOsl 



dt 



Vu) = diffusion + emission — absorption, 



(5) 



which contains the plasma velocity u of the background plasma. Here the colon denotes the contraction 
of the two tensors Pj, and Vu. The processes described by the symbolic terms on the right hand side of 
equation ([S]) will be described below. 
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Setting the Eddington factor = 1/3 corresponds to the radiation diffusion model. In this case the 
radiation is assumed to be effectively isotropic and the spectral radiation pressure can be described by the 
scalar 

P. = ^E, = hr-l)E,, (6) 

where we have introduced the adiabatic index of the radiation field, which in this case has the relativistic 
value 7r — 4/3. The time evolution for the spectral energy density can then be simplified to 

dE dE 

+ ^ ■ (^fU) — — 1)(V • \i)v ^ ^ = diffusion + emission — absorption. (7) 

The second and third terms on the left hand side of equation ([7]) express the change in the spectral energy 
density due to the advection and compression of the background plasma, which moves with the velocity u, 
as well as the frequency shift due to compression. In the free-streaming limit where the radiation hardly 
interacts with the matter, will approach one. In this paper we will keep Xi^ = 1/3 Sind at the same time 
use a flux limited diffusion for the free-streaming regime whenever needed (see Section 12. 5[) . 

The set of equations for the spectral energy density ([7]) still consists of an infinite amount of equations, 
one for each frequency. A finite set of governing equations to describe the radiation transport in the multi- 
group diffusion approximation is obtained when we choose a set of frequency groups. Here we enumerate 
groups with the index, g — 1,...,G. The interval of the photon frequencies, relating to the f;th group 
is denoted as [1^9-1/2, '^9+1/2]- A discrete set of group energy densities. Eg, is introduced in terms of the 
integrals of the spectral energy density of the frequency group interval: 

Eg= E^dv. (8) 

Now we can integrate equation Q to arrive at the desired set of the multigroup equations: 

+ ^ • (^9") + (> - 1)^5^ • u - (7, - i)(v • u) r^^'^ ^-^^dv 



= / (diffusion -I- emission — absorption)^!/. (9) 

The fourth term on the left hand side is a frequency shift due to the plasma compression. This term is 
essentially a conservative advection along the frequency axis. 

In the context of the multi-group radiation diffusion, a discussion about the stimulated emission is not 
less important than LTE. E xcellent accounts on the stimulated emission exist in the literature, see for instance 



Zel'dovich fc Raizer I ()2002l) . Here, we merely sumarize how the stimulated emission modifies the absorption 
opacity obtained from, e.g., opacity tables. This is important when dealing with externally supplied 
opacity tables, since the CRASH code assumes that the absorption opacities are corrected. Integrating the 
total absorption and emission over all directions and summing up the two polarizations of the photons, the 
following expression can be derived for the emission and absorption 

emission — absorption = ck°' {Bi, — E,y) , (10) 

where the effective absorption coefficient, is introduced to account for the correction due to stimulated 
emission: 

/ r 1 \ 

(11) 



/c^' = I 1 - exp 
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in which e = is the photon energy, ks is the Bohzmann constant, and Tg is the electron temperature. 
We also introduced the spectral energy density distribution of the black body radiation (the Planckian) 



Stt 



(12) 



The total energy density in the Planck spectrum equals B ~ dvBi, = aT^, where a = Sn^ kg/ {15 h^c^) is 
the radiation constant. 

We use the standard defini tion of the group Planck mean opacity Kpg and group Rosseland mean opacity 
KRg iMihalas fc Mihalas lEosi ) 



dBg 



dvBy 



(13) 



Vg-l/-2 



in which k*, is the spectral total op a city. The righ t hand side of equation ([9]) can now be written as (see for 



instance 



Mihalas fc Mihalas I (|l984l ): IPomraming I (|2005f l) 



dt 



V • [EgVL) + (7, - l)EgV • u - (7, - 1)V • u /" 

•^I'g-l 

V ■{DgVEg)+ag{Bg~Eg), 



^3 + 1/2 d{vE^) 



dv 



(14) 



where Dg = c/{3Kug) is the radiation diffusion coefhcient for radiation group g. The absorption and emission 
uses the coefficient Cg — cnpg. These group mean opacities are either supplied by lookup tables or by an 
opacity solver. 

In a single group approximation (gray diffusion) , the spectral energy density is integrated over all photon 
frequencies and the total radiation energy density is obtained by 



Er{v,t) 



Eudv. 



(15) 



This amounts to su mming up all groups Er = Y'. g Eg. The g r ay ra di ation d i ffusion equation can be derived 
as (see for instance iMihalas fc Mihalas I (|l984h : IPomraming I (|2005l ): [Prakel ^20o4 )) 



dEr 



V • {ErVL) + (7r - l)ErV • U = V • {DrVEr) + CTr(S - Er 



(16) 



where the diffusion coefficient is now defined by the single group Rosseland mean opacity kji as D,. — 
c/ {itipi), and the absorption coefficent ar is defined in terms of the single group Planck mean opacity Kp as 

(Jr — CKp. 



2.2. Hydrodynamics 

In the CRASH code, a single fluid description is used, so that all of the atomic and ionic species as well 
as the electrons move with the same bulk velocity u. The conservation of mass 

|^+V.(pu)=0, (17) 
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provides the time evolution of the mass density p of all the materials in the simulation. The plasma velocity 
is obtained from the conservation of momentum 

^ + V-[puu + /(p + p,)] = 0. (18) 

The total plasma pressure is the sum of the ion and electron pressures: p = pi + Pe- The net force of 
the radiation on the plasma is given by the gradient of the total radiation pressure —Vpr, where the total 
radiation pressure is obtained from the group radiation energies: pr = {'jr — 1) X^^s- 

In a high density plasma, the electrons are very strongly coupled to the ions by collisions. However, 
for higher temperatures, the electrons and ions get increasingly decoupled. At a shock front, where ions are 
preferentially heated by the shock wave, the electrons and ions are no longer in temperature equilibrium. 
Ion energy is transfered to the electrons by collisions, while the electrons will in turn radiate energy. We 
therefore solve separate equations for the ion/atomic internal energy density Ei and the electron internal 
energy density i?e: 



dt 
BE, 
dt 



V-(£;,u)+p,V-U=f7,e(re~T,), (19) 

G 

V • {E.vl) +peV • u = V • (CeVTe) + C7,,{T, - T,) ag{Eg - Bg). (20) 

3=1 



The coupling coefficient aie — riakB/Tie in the collisional energy exchange between the electrons and ions 
depends on the relaxation time Tie{Te, Ua, m) and the atomic number density Ua- The energy transfer depends 
also on the difference between the ion temperature Ti and the electron temperature Tg. In equation (|20p . 
we have included the electron thermal heat conduction with conductivity Ce(Te, n^, m). Since the electrons 
are the species that are responsible for the radiation absorption and emission, the energy exchange between 
the electrons and the radiation groups is added to equation (PO]) . 

For the development of the numerical schemes in Section [3] we will use an equation for the conservation 
of the total energy density 

2 G 
e = ^+E,+E,+Y,Eg, (21) 

9=1 

instead of the equation for the ion internal energy ([T9|) . This is especially important in regions of the 
computational domain where hydrodynamic shocks can occur, so that we can recover the correct jump 
conditions. The conservation of the total energy can be derived from equations (fT4| and ([T7| - ((20|) as 

^ + V • [(e + p + pr)n] = V • (CeVTe) + ^ V • {DgVEg) . (22) 

9=1 

The frequency shift term in equation (jl4p due to the plasma compression does not show up in the conservation 
of the total energy if we use energy conserving boundary conditions at the end points of the frequency domain, 
i.e. at = and v ~ oo m the analytical description or at the end points of the numerically truncated finite 
domain. 



2.3. Level Sets and Material Identification 



In many of the CRASH applications, we need a procedure to distinguish between different materials. 
We assume that the materials do not mix, but differ from each other by their properties such as the equation 
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of state and opacities. If we use M different materials, tlien we can define for each material m = 1, . . . , M 
the level set function dm{'r,t) that is initially set to zero at the material interface, while positive inside the 
material region and negative outside. Generally, we use a smooth and signed distance function in the intial 
state. At later times, the location of material m is obtained by means of a simple advection equation 

dd 

^ + V • (rf„u) = d,„V • u. (23) 

For any given point in space and time, we can determine what the material is, since analytically only one 
of the level set functions dm can be positive at any given point. Numerical errors will create regions where 
this is not true. In practice we take the largest dm- This is a simple approximation, we may explore more 
sophisticated approaches in the future. The number of material levels M can be configured at compile time. 



2.4. Equation of State and Opacities 

We have implemented EOS solvers and a code to calculate the frequency averaged group opacities. The 
implemention will be reported elsewhere, but mention that in the EOS and opacity solver the temperature 
is assumed to be well below the relativistic values: T <C 10^ eV. The non-relativistic speed of motion is also 
assumed while simplifying the radiation transport equation and to neglect the relativistic corrections for 
opacities. In this paper, we will assume that all necessary quantities are calculated and stored in lookup 
tables. Our EOS solver assumes that the corrections associated with ionization, excitation, and Coulomb 
interactions of the partially ionized ion-electron plasma are all added to the energy of the electron gas and 
to the electron pressure. This is possible since those corrections are controlled by the electron temperature. 
The ion internal energy density, ion pressure, and ion specific heat in an isochoric process per unit of volume 
are simply 

Pi = riakBT,, Ei = -, Cvi = = -, (24) 

7-1 p 7-1 

which are due to the contributions due to ion translational motions for which 7 = 5/3. 

The relations among the electron internal energy density, pressure, density, and temperature are known 
as the EOS. To solve these relations is usually complex and time consuming. We therefore store these relations 
in invertible lookup tables. For each material m, our EOS tables have the logarithmic lookup arguments 
(logTe, log ria). The list of quantities stored in these tables is indicated in Table [TJ These lookup tables are 
populated with quantities that are needed for both single temperature and two temperature simulations. 
For two-temperature plasma simulations, we will need Pe, E^,, the electron specific heat Cve, the electron 
speed of sound gamma 75^. For convenience we added the total pressure p ^ + Pi, total internal energy 
density E = Ei+ E^, single temperature specific heat Cy, and the single temperature speed of sound gamma 
75, which can be used in single temperature simulations. We use high enough table resolutions so that it 
is sufficient to use a bilinear interpolation in the lookup arguments. If Pe or E,, (or p and E in single 
temperature mode) are known on entry of the lookup instead of Te, we do a binary search in the table to 
find the appropriate electron temperature. The latter only works as long as the necessary thermodynamic 
derivatives are sign definite, i.e. the table is invertible. Other thermodynamic quantities that are needed, 
but not stored in these lookup tables, can be derived. For example, the electron density can be obtained 
from the mean ionization rig = UaZ. 

In addition, we have lookup tables for the averaged multigroup opacities. These tables are either 
constructed internally for a given frequency range, number of groups, and the selected materials, or externally 
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supplied. For any material m, the logarithmic lookup arguments are (logp, logTg). The stored quantities, 
see Table [21 are the specific Rosseland mean opacity Kjig/ p and the specific Planck mean opacity Kpgj p for 
all groups g — 1 , . . . , G that are used during a simulation. The Planck opacities are assumed to be corrected 
for the stimulated emission, as discussed in Section [2. II The groups are always assumed to be logarithmically 
distributed in the frequency space. 



2.5. Flux Limited Diffusion 



In 



Radiation diffusion theory can transport energy too fast in the optically thin free streaming limit 
the diffusion limit, the radiation diffusion flux for each group follows Pick's law Fg = —DgVEg, where the 
diffusion coefficient Dg depends on the Rosseland mean opacity Kjig for the group g via Dg = c/{3Kfjg). 
This flux is however not bounded. In the optically thin free-streaming limit, the magnitude of the radiation 
flux can be at mos t cE„ i n order to mainta i n cau s ality. Various flux limiters exist in the literature, see for 



instance 



Minerbo I (|l978l ); iLund fc Wilson I (|198Clf ): iLevermore fc Pomraming I (|l98lh . that ensure that the 



diffusion flux is limited by this free streaming flux. We implement ed the so-cal led square-root flux limiter 
to obtain the correct progation speed in the optically thin regime dMorellboOoh . For this flux limiter, the 
diffusion coefficient is rewritten as 

Dg = ' (25) 



In the limit that the radiation length scale Ln = Eg/\WEg\ is large, the diffusive limit is recovered. For a 
small radiation length scale, Dg = c\Eg\/\W Eg\ and the radiation propagates with the speed of light. 



Similarly, we implemented the option to limit the electron thermal heat fiux (see Drake ( 2006h for more 
details on electron flux limiters). The classical Spitzer-Harm formula for the colhsional electron conductivity 
is proportional to Te^^/Z'^, where is the mean square ionization of the used material. The colhsional 
model is only valid when the temperature scale length Lt = Te/|VTe| is much larger than the colhsional 
mean free path of the electrons Xmfp- When the temperature scale length is only a few Xmfp or smaller, 
this description breaks down. This may for instance happen in laser-irradiated plasmas. One could in that 
case find the heat fiux from solving the Fokker-Planck equation for the electrons, but this is computationally 
expensive. Instead, we use a simplified model to limit the electron heat flux. A free-streaming heat flux can 
be defined as the thermal energy density in the plasma transported at some characteristic thermal velocity: 
Fps = neksTeVth, where vth = \fk^T^Jjn'e. For practical applications, the maximum heat transport is 
usually only a fraction of this free-streaming flux: F = — (/i^i?5/|VTe|)VTe, where / is the so-called fiux 
limiter. T his heat fiux model is the threshold model and is also used in other radhydro packages, such as 
HYADES (jLarsen &: Lane 1119941 1. The flux- limited heat fiux can now be defined as 



F^-min(C.,^)vT. 



(26) 



The flux limiter / is an input parameter and can be tuned to let the simulated results better fit reality. 



3. THE NUMERICAL METHOD 



In this section, we describe the discretization of the set of multi-material, radiation hydrodynamics 
equations for the density pT)) . momentum (fTS)) . total energy ((22)) . electron internal energies ((20)) . radiation 
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group energy ([M]), and material level set functions (|23l) . The equations are time integrated using an operator 
split method to solve the equations in substeps. Formally, we may write this system as 

-qI = %ydro(U) + R'frequency(U) + R-diffusion(U), (27) 

where U is the vector of state variables. We have split the right hand sides of the equations into three parts 
and time advance the equations with an operator splitting method in the following order; (1) The right hand 
side Rj^yjj,Q describes the advection and pressure contributions (Section 13. ip . This part is essentially the 
ideal hydrodynamic equations augmented with the advection and compression of the radiation energy, the 
electron internal energy, and the material level sets. (2) The right hand side R-frequency advection 
of the radiation field in frequency space (Section 13. 2|) . (3) The right hand side R-ciiffusion t^^kes care of 
the diffusion and energy exchange terms, which we will solve with an implicit scheme (Section 13.31) . The 
operator splitting is not unique. Instead of splitting the hydrodynamic advection operator and the extra 
advance operator for the frequency advection, one may attempt to discretize the frequency advection flux 
as an extra flux for the control volume of the four-dimensional (x, y, z, v) space. However, since the CRASH 
code is built around the existing BATS-R-US code in ID, 2D, and 3D, we opted for splitting the frequency 
advection from the hydro update. The boundary conditions are treated in Section [2111 



3.1. Hydro Solve 

In the first step of the operator splitting, we update the hydrodynamic equations including the advection 
and compression of the radiation energy density, electron internal energy density and the level sets. We have 
implemented two variants to solve the hydrodynamic equations: (1) using conservation of the total energy 
(Section I3.1.ip and (2) a non-conservative pressure formulation (Section I3.1.2p . We can combine the two 
discretizations in a hybrid manner a simulation. 



3.1.1. Conservative 



We have implemented several hydrodynamic shock-capturing sche mes in the C RASH code: the HLLE 



scheme (Harten et al. lll983l : lEinfeldt et al. Ill99ll) . the Rusanov scheme (|Yee 1119891 ) . and a Godunov scheme 



dCodunov II 1959!) ^ with an exact Riemann solver. In this section, we will explain how we generalized the HLLE 
scheme for our system of equations that includes radiation, level sets, and an EOS. The other hydrodynamic 
schemes can be generalized in a similar fashion. 

Typical hydrodynamic solvers in the literature assume constant 7. Our problem is to generalize the 
constant 7 hydro solvers for the case of spatially varying polytropic index, 7e, due to ionization, excitation 
and Coulomb interactions. A method that is applicable to all the aforementioned, constant 7, hydrodynamic 
shock-capturing schemes is to split the electron internal energy density as the sum of an ideal (transla- 
tional) energy part Pe/(n ~ 1) and an extra internal energy density Ex- Similarly, we can define an ideal 
total energy density 

ei 



pu 



Pi +Pe 

7-1 



G 



(28) 



which is related to the total energy density hy e = ej + Ex ■ We will time advance Pe with the ideal electron 
pressure equation and Ex by a conservative advection equation, and then apply a correction step as described 
below. 
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The time update with the operator Rj^yj^j-Q solves the following equations: 

V • ipu) = 0, (29) 



dp 
dt 



^ + V-[puu + /(p + p,)]=0, (30) 

dei 
~dt 



V- [(e/+p + p,)u] =0, (31) 



^ ' ^ V-(peU)+PeV-U = 0, (32) 



7 — 1 9t 7 — 1 

dEx 



dt 

dEg 

dt 
ddm 
~dt 



V • [Exn] = 0, (33) 

V • [Egu) + (7, - \)EgV • u = 0, (34) 

V • (d„u) - d,„V • u = 0, (35) 



where the frequency advection, diffusion, and energy exchange terms are ommitted in this first operator 
step. After each time advance from time to time we have to correct e, e/, pe, and Ex- We denote 

the uncorrected variables with a superscript *, then we recover at time level n + 1 the true electron internal 
energy E'"'*'^ and the true total energy density e"''^^ by 

^„+i ^ _P^^E*x. (36) 
7- 1 

e"+i = e*i + E*x. (37) 

Since both e/ and -Ejf follow a conservation law, the total energy density e is also conserved. The true 
electron pressure is recovered from the updated electron internal energy and mass density by means of the 
EOS: 

pr'=PEos(p"+\i?r\™), (38) 

where the function Peos can be either a calculated EOS or an EOS lookup table for material to, determined 
by the level set functions rf^"^^ fSection l2.3p . The extra internal energy Ex is reset as the difference between 
the true electron internal energy and the ideal electron internal energy for 7 = 5/3: 

S^+i = E'^+^ - ^ — . (39) 
7-1 

This is postive because the EOS state peos satisfies E,. — Pe/il — 1) > at all times. The ideal part of the 
total energy density at time level n + 1 can now be updated as 

^g„+i (40) 

We have now recovered e""*"^, e/""*"^, p""''^; ^^^d E^"^ at time 

We time advanced the hydrodynamic equations to the time level * with a shock-capturing scheme with 
a constant 7 — 5/3. For an ideal EOS, the speed of sound of the equations (|29l) - (p4|) can be derived as 



hjPi + Pe) + IrPr ^^^-^ 
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which includes the modifications due to the presence of the total radiation pressure. This speed of sound will 
be used in the hydro scheme below. Since the CRASH EOS solver always satisfies Ex > and 7e < 5/3, 
the speed of sound for the ideal EOS is always an upper bound for the true speed of sound. 



We use shock-capturing schemes to advance the equations (129)) -([35 |) . In the following, we denote the 
(near) conservative variables hy U — {p, pu, ej,pe, Ex, Eg, dm) and let U be grid cell averages in the standard 
finite volume sense. If we assume for the moment a ID grid with spacing Ax, cell center index i and cell 
face between cell i and i + 1 identified by half indices i + 1/2, then we can write the two-stage Runge-Kutta 
hydro update as 



1+1/2 



2Ax (•^^+1/2 "^'"-1/2) 



n+l 



2Ax 
At 



n+1/2 j.n+1/2 
J 7 



-1/2 



(42) 
(43) 



where / is the numerical fiux. In particular, the HLLE fiux / equals the physical fiux ^(t^/^1/2) when 



Ui + Cs < 0, i^(?/i^i/2) ^^en 



Cs > 0, and in all other cases it uses the weighted flux 



/i+l/2 

Here, the left and right cell face states are 



-(ttR 

i+l/2 



u 



-1/2) 



ct 



Cs 



1/2 

TjR 

'^i+1/2 



We use the generalized Koren limiter, and define the limited slopes as 



A Ui — minmod 
A^Uj = minmod 



/3(C/,+i -[/,), /3(C/.-C/.-i), 



/3(C/,+i -f/,),m-C/.-i). 



(44) 

(45) 
(46) 

(47) 
(48) 



for the extrapolatio ns from the l eft and right. Th is reconstruction can be third order in smooth regions 
away from extrema ( Koren 19931 : Toth et al. 2008 ). The parameter j3 can be changed between 1 and 2, but 
in simulations with adaptive mesh refinement we have best experience with /3 = 3/2. We generally apply 
the slope limitcrs on the primitive variables (p, u,pi,pe, Ex/ p, Eg, dm), instead of the conservative variables. 
We apply the slope limiter on Ex/ p instead of Ex since Ex/ p is smoother at shocks and across material 
interfaces. A multi-dimensional update is obtained by adding the fluxes for each direction in a dimensionally 
unsplit manner. 

After each stage of the two step Runge-Kutta, we correct for the EOS effects via the update procedure 
outlined in equations ([Mt - ipO]) . 



3.1.2. Non-Conservative Pressure Equations 



In regions away from shocks it is sometimes more important to preserve pressure balance than to have a 
shock capturing scheme that recovers the correct jump conditions. This is especially important at material 
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interfaces. We therefore have implemented the option to solve the hydro part of the pressure equations 

^ + V-(p,u) + (7-l)p,V-u = 0, (49) 

^ + V • (peU) + (7Se - l)PeV • U = 0, (50) 

instead of the equations for the total energy (pij) and the electron internal energy (j3D) . As long as the speed 
of sound gamma for the electrons 

IS. = ^ f (51) 

Pe\dp J 

is smaller than 7 = 5/3, the numerical scheme is stable. Contrary to the energy conserving scheme, the 
pressure based scheme can directly include the EOS and we no longer need the time evolution of the extra 
internal energy density p3|) . The EOS contribution in the electron pressure equation (|50|) is implemented 
as a source term —{'fse — 7)PeV • u added to the ideal electron pressure equation. 

To facilitate using both the shock capturing properties and the pressure balance at the material interfaces 
during CRASH simulations, we have several criteria to automatically switch between them accordingly. One 
of the criteria, for instance, uses a detection of steep pressure gradients as a shock identification. The user 
can select the magnitude of the pressure gradient above which the scheme switches to the conservative energy 
equations. 



3.2. Frequency Advection 

The set of multigroup equations (|14p contain an integral over the group photon frequencies. Performing 
this integration, the frequency advection update by the R-frequency operator can be written as 

dE 

- {jr ~ 1)(V • U) [l^g+l/2E,,{iyg+i/2) - l^g-l/2E^{l^g-i/2)] = 0. (52) 

These equations, however, do still depend on the unknown photon frequency v and the spectral radiation 
energy density E^. We will now restrict the analysis to a frequency grid that is uniformly spaced in the 
frequency logarithm, i.e., 

ln(z/g^i/2) ~ lii(t'g_i/2) — A(lni/) = constant. (53) 

For large enough number of frequency groups G, the group energy Eg can then be approximated as the 
product of the photon frequency, spectral radiation energy E^^ and the logarithmic group spacing A(lni/): 

Eg^ I E^dv^ i E^vd{\nv) Eyvl^{\nv). (54) 

Using this approximation in equation (j52p . we obtain our final form of the frequency advection 

dEg Eg+l/2 - Eg_l/2 



dt " A(lnz/) 



0, (55) 



where u,, — — (7r — 1)V • u is the frequency advection speed. The values i?g±i/2 are to be interpolated from 
the mesh-centered values Eg towards the group boundaries. 

The frequency advection is a conservative linear advection in the log-frequency coordinate, for which 
the physical flux is defined as Fg_i/2 = u^Eg_i^2- For the boundary conditions in the frequency domain we 
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assume zero radiation flux so that no radiation will leak at the edges of the frequency domain. Equation 
(|55p can be discretized with the one-stage second order upwind scheme 



fg+1/2 - fg-1/2 



A (In:/) 

where time level * is now the state after the hydro update and the numerical flux is 



/g-l/2 
/g-l/2 



1-C 



-A{Eg+l - Eg, Eg " Eg^l] 



En - 

^ 2 



1 — C - 

Eg-1 H A{Eg - Eg^l,Eg^l ^ Eg^2) 



Ul> < 0, 

> 0. 



(56) 



(57) 



and we use the superbee limiter ( Roe 19861) for the limited slope A. The Courant-Friedrichs-Levi number 
C — Ai/A(ln ^) depends on the hydrodynamic time-step At. If C is larger than one, the frequency 
advection is sub-cycled with the number of steps equal to the smallest integer value larger than C. 



3.3. Implicit Diffusion and Energy Exchange 

The stiff parts of the radhyro equations are solved implicitly in an operator split fashion. These stiff 
parts are the radiation energy diffusion, electron heat conduction, and the energy exchange between the 
electrons and each energy group g and between the electrons and ions. In this section, we will describe two 
implicit schemes that are implemented: (1) solving for all radiation groups, electron and ion temperatures in 
a coupled manner (Section l3.3.1|) . and (2) solving each radiation group energy and the electron temperature 
independently (Section l3.3.2p . Our strategy for resolution changes is described in Appendix El while the 
modifications for the rz-geometry are explained in Appendix [Bl 



3.3.1. Coupled Implicit Scheme 

Discretizing the diffusion and energy exchange terms of equations ([T9| -(|20 l) . and ([T4l) implicitly in time 
leads to 

<(T;+i-Tf+i), (58) 

G 

al{T^+^ - Tr+i) -I- V • C:VT:+^ + ^ al{E^+^ - B^^+^), (59) 

9=1 

a;(i?;^+i - i?;'+i) + v • dive'^^+\ (60) 

where time level * now corresponds to the state after the hydro update and the frequency advection. The 
coupling coefficients cr*, and cr* and the diffusion coefficients C* and D* are taken at time level * (frozen 
coefficients). One can either (1) solve the coupled system of G 4- 2 equations ([58|) -(|60 l) implicitly or (2) solve 
equation ([58]) for the ion internal energy E"^'^ , substitute the solution back into equation ((59|). and solve the 
resulting reduced set of G -I- 1 equations ((5^ - (l5(Il) implicitly. Here we describe the second scheme, because 
it is more efficient, especially for small number of groups, e.g., for gray radiation diffusion. Note that if we 
did include ion heat conduction in (1601) . then we would have to solve the entire coupled system of equations. 



E^+^ - E* 



At 




En+l _ 


e; 


At 






e; 
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First we introduce the ion Planck function Bi = aT^ as a new variable similar to the electron Planck 
function B = aT^, and replace Ei and with these variables using the chain rule 



dt 



dE, dT, dB, Cv^ dB, 



dT, dB, dt AaTf dt 



dE, 
dt 



Cve dB 
Aan~dt' 



(61) 



in which Cvi and Cve are the specific heats of the ions and electrons, respectively. Now equation (|58|) can 
be replaced with 

B-+' = B* + AiaL(S"+^ - B-+'), (62) 



where 



4aT3 



1 



Cv^ a{T, + T,){T^+Tl)^ 



(63) 



is again taken at time level *. The numerator comes from (T^ — T^)/{Te — Ti). Equation (p^ can be solved 
for This result can be substituted into the electron internal energy equation ((59)) to obtain 



At 



_ B*) = (j'l^{B* - B"+i) + V • C^VS"+i + (r*g{E^^^ - w;B"+1), 



where we have introduced new coefficients at time level 



Vi 



AaTf 1 + Ma'„ 



C' = 



c: 



AaTr 



(64) 



(65) 



The Planck weights — B*/B* satisfy J2g''^g — 1- It is convenient to introduce the changes AB 



At 



a" - V • cy 



AB 



J2 S*(Ai?3 - w;AB) = aUB* - B* 

3=1 

G 

^■cyB* + j2'^;{E;-w;B*), 



— - V • D*,V 
At » 



AE„ 



9=1 

<j:{w:AB - AEg) = aMB* - E*) + V • DyE* 



(66) 
(67) 



This is a coupled system of G + 1 linearized equations for the changes AB and AEg . The right hand sides 
are all at time level *. 

A discrete set of equations are obtained by applying the standard finite volume method to the equations 
()66p and ()67|) and partitioning the domain in a set of control volumes Vi, enumerated by a single index 
i — 1, . . . , I . As an example, the fluxes Fgij associated with the radiation diffusion operator may be obtained 
by approximating the gradient of the group energy density with a simple central difference in the uniform 
part of the mesh: 

- / V • iDgVEg)dV = J2 ^a^^ = E ^'^^g'^ ^'lf'i ' (68) 

where the index j enumerates the control volumes which have a common face with the control volume i, 
the face area being Sij, and the distance between the cell centers is |xi — Xj|. N ote that we assume d here 
an orthogonal mesh. Generalization to curvilinear grids can be done as shown in lToth et al. ( 2008 ). The 
diffusion coefficients at the face are obtained by simple averaging of the cell centered diffusion coefhcient: 
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Appendix [X] 



Dgj)/2. The discretization of the diffusion operator at resolution changes is described in 



The hnear system (p^ - dST)) can be written in a more compact form as the hnearized imphcit backward 
Euler scheme 

'l-Ai|^^ AU = A<R(U*), (69) 

where U are the / x (G + 1) state variables B and Eg for all / control volumes, and AU = U"+^ — U*. 
R is defined by the spatially discretized version of the right hand side of equations and (p7|) . The 
matrix A = I — AtdR/dXJ is a / x / block matrix consisting of (G + 1) x (G + 1) sub-matrices. This 
matrix A is in general non-symmetric due to the Planck weight w* in the energy exchange between the 
radiation and ele ctrons. To solve this s ystem we use Krylov sub-space type iterative solvers, like GMRES 
(|Saad fc Schultz Ill9 86) or Bi-CGSTAB (van der Vorst II1992I) . To accelerate the convergence of the iterative 
scheme, we use a preconditioner. In the current implementation of CRASH, we use the Block Incomplete 
Lower-Upper decomposition (BILU) preconditioner, which is applied for each adaptive mesh refinement block 
independently. For gray radiation diffusion the Planck weight is one and the matrix A can be proven to 



be sym metric positive definite (SPD) for commonly used boundary conditions (see for example [Edwards 
( 1996 )'). I n tha t case we can use a preconditioned conjugate gradient (PCG) scheme (see for instance 



Eisenstat I (|198l[) ). 



For some verification tests, we can attempt to go second order in time under the assumption of temporally 
constant coefficients using the Crank-Nicolson scheme 



U 



n+l 



_ u* 



At 



(l-a)R(U*) + aR(U 



n+l\ 



(70) 

(aR/9U)*AU to obtain 
(71) 



with a = 1/2. The implicit residual can again be linearized R(U"+^) = R(U* 
the linear system of equations 

I - aAt|^^ AU = AtR(U*). 

We use the same iterative solvers as for the backward Euler scheme. 

Finally, we show how we use the solution AB and AEg for g = 1,...,G from the non-conservative 
equations ([55)1 and (p7|) to advance the solution of the original equations ([551) - (pn)) and still conserve the 
total energy. One needs to express the fiuxes and energies on the right hand side in the latter equations in 
terms of and E"^^ while still keeping the coefficients frozen. After some algebra we obtain 



- E* + AtaUB"^+' - B*), 

= £;:+Gt.,(B"+i-B*), 

= E*+AEa. 



(72) 
(73) 
(74) 



This update conserves the total energy to round-off error. Note that at this final stage, taking too large time 
step may result in negative ion internal energy E"~^^ if ^ B* and negative electron internal energy 

^n-i-i ^n-i-i <^ B* . If this happens, the advance might be redone with a smaller time step, to limit the 
drop in B, or by some other timestep control scheme. A generalization of the conservative update to the 
Crank-Nicolson scheme is also implemented for verification tests with time constant coefficients. 

For completeness, we mention that in the absence of radiation we solve during the implicit step for 
the temperatures and Ti instead of the radiation energy-like variables aT^ and aT,^. In that case the 
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corresponding matrix A is always SPD. In princi ple, the formulation in temperatures can be generalized 



to radiation as well. In Landau fc Lifshitz ( 1980f) . a spectral temperature T^{E^,v) is defined, such that 



the spectral energy density is locally equal to the spectral Planckian energy density at the temperature Ty: 
Eu = B^{T^, v). This relationship is a one-to-one map. A group temperature, Tg, can also be introduced as 
the discrete analog such that the group energy density can be obtained by 



B^{Tg,y)dv. 



(75) 



The equation (|60|) can be recast as equation for the group temperature Tg. This introduces the group specific 
heat of the radiation Cg = dEg/dTg. The set of equations (|58 l) - (|60l) reform ulated as an imp licit backward 



Euler scheme for the temperatures Ti, Te, and Tg can in a similar way as in lEdwards I (|l996[ ) be proven to 
be SPD. While this scheme has the advantage of being SPD, the conservative update of the group energy 
density E^^^ = E* + C*ATg might result in negative energy density E"^^ for too large time steps. 



3.3.2. Decoupled Implicit Scheme 



The coupled implicit scheme of Section [3.3.11 requires the solution of a large system of equations [G + 1 
variables per mesh cell). The preconditioning for such a system can be computationally expensive and 
requires overall lots of memory. We therefore also implemented a decoupled implicit scheme that solves each 
equation independently. 

For some applications, the electron temperature does not change much due to energy exchange with the 
radiation. This is typically so if the electrons have a much larger energy density than the radiation, so that 
Tf. changes little due to the interaction with the radiation in a single time step. In that case, we solve first for 
the electron and ion temperatures without the contributions from the radiation-electron energy exchange. 
Let again time level * indicate the state after the hydro update and frequency advection, and freeze again 

c: 



a* at time level *. Discretization in time now leads to 



El 



■n+1 



E* 



At 



At 



t: 



'n+l\ 



V • c:vt: 



(76) 
(77) 



where the time level ** of E'e indicates that we still have to do an extra update to time level n + 1 with the 
radiation-electron energy exchange. Each radiation group energy density is solved independently using time 
level * for the electron temperature in B* : 



-'^aiBg-Eg ) + 

where we have exploited the assumption that B* is not stiff. 
Equations ((75)) - ([75)) can be recast in equations for the G - 

AEg ^ E^+^ - e;-. 



d;ve-+\ 



(78) 



1 independent changes AB = B** - B* and 



At 
1 

At ' 



a'L - V • C:V 



AB = a'l^{B* -B*) + V ■ G'^VB* 



AEg = 



a;iw;B*-E;) 



W -Dl^E*. 



(79) 
(80) 
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where we have used the definitions (|6ip . ((63)) . and (1651) of the coefficients, frozen at time level *. Each 
equation for the changes is in the form of the linearized implicit backward Euler scheme (j69l) and can be 
solved independently with iterative solvers like GMRES and Bi-CGSTAB using a BILU preconditioner. 
As long as the boundary conditions are such that the matrices are symmetric and positive definite, the 
preconditioned conjugate gradient method might also be used. 

In a similar manner as with the coupled implicit scheme, a conservative update for the energy densities 
can be derived as 

E^+' = E;+AEg, (81) 
= E* +AtaUB** - B*), (82) 

G 

£;«+! = E*^+C{,^{B**^B*) + AtY,'^*g{Eg+^-w*gB*), (83) 

9=1 

that preserve the total energy to roud-ofF errors. The main difference between the conservative update in 
the coupled and decoupled schemes is that here the energy exchange between the radiation and electrons is 
added afterwards as the last term in equation (1551) . 

This scheme requires less computational time for preconditioning and the Krylov solver than the coupled 
implicit algorithm, however it generally needs more message passing in parallel computations. It is therefore 
not always guaranteed that the decoupled scheme is faster. The memory usage is always smaller. 



3.4. Boundary Conditions 

The CRASH code allows for any user specified type of boundary conditions. Several commonly used 
boundary conditions are readily available in the main code for convenience, e.g., fixed, extrapolation with 
zero gradient, periodic, and reflective boundary conditions. 

For the radiation field, we have implemented a zero or fixed incoming fiux boundary condition that is 
used instead of the extrapolation with zero gradient. This type of boundary condition is useful if there are 
no sources of radiation outside the computational domain and we assume that outflowing radiation does not 
return back into the computational domain (zero albedo). Note that simple extrapolation with zero gradient 
can make the radiation diffusion problem ill-posed. The boundary condition is derived as follows: Radiation 
diffusion approximation corresponds to a linear-in-angle intensity distribution 

Ig - ^E, + Af, . n, (84) 

so we can calculate the radiation flux through a boundary surface. If we define the outward pointing normal 
vector of the boundary as n^, the net fiux of radiation energy inward through this boundary is 

= - / n,- nigdn = ^ - in, . F„ (85) 

where the closure (|M)) is used. In the radiation diffusion model, the flux is written as Fg — —DgWEg, where 
the diffusion coefficient Dg is a nonlinear function of Eg and V Eg in a flux limited diffusion model. The 
boundary condition satisfies 



2Dn 4 ^ 

9 



Eg + ^nb-VEg = -Fl". (86) 
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For the left boundary in the x-direction, for instance, this can be discretized as 

EgO + Egl 2Dg Egl — E gQ _ 4 



f;", (87) 



where the index 1 corresponds to the last physical cell and to the ghost cell. This equation can be solved 
for the ghost cell value. For zero incoming radiation flux boundary conditions we set F™ = 0. 



4. CODE VERIFICATION 

To test the CRASH as well as the BATS-R-US and SWMF codes, we have implemented numerous tests. 
These tests are subdivided in two categories: functionality tests and verification tests. Both test suites are 
performed automatically and return pass or fail messages depending on whether or not certain predefined 
tolerance criteria are met. This automated testing process provides a software quality confidence especially 
when used in combination with a software version control system like CVS (Concurrent Versions System) to 
recover previous correctly performing code. 

The functionality tests are performed nightly on several computer platforms with different compilers 
and number of processors. They consist of unit tests and full system tests. The unit tests are designed to 
test a particular unit, for example a linear equation solver. The full system tests on the other hand, exercise 
the code in the way the end-users will use it for their research applications. We always try to cover as much 
code as possible with these tests so that we can discover bugs and other unwanted side effects early on. 

To test the correctness of the implemented algorithms we have also constructed a suite of verification 
tests. This suite is executed daily on a dedicated parallel computer and runs specific simulations to quantify 
against analytic and semi-analytic solutions, whenever possible. The CRASH test repository currently covers 
a wide range of tests for hydrodynamics, multi-material advection methods, gray and multigroup radiation 
diffusion, heat conduction, to mention a few. These are performed to test for grid and/or time convergence, as 
deemed necessary. We also simulate full system laboratory experiment configurations in various geometries, 
dimensionality, and physics fidelity. The results are either validated against laboratory experiments or simply 
used to check that the code keeps performing these simulations as expected. Once a week, we also perform a 
parallel scalability test on a large parallel computer to verify that the code does not degrade in performance 
during further development of the software. 

In the following sub-sections, we highlight some specific verification tests related to the implicit radiation 
(Section l4.2l) and heat conduction fSection l4.3p solver. The tests cover both Cartesian and rz-geometry, and 
some of them also involve the hydrodynamic solver. We demonstrate a 3D full system test in Section 14.41 
and describe the parallel scalability in Section H3] 



4.1. Error Assessment 

For the assessment of the accuracy of the solutions in the test suites, an appropriate definition of the 
numerical errors have to be defined. We will use two types of errors to quantify the verification analysis: 
The relative LI error is defined as 

a=l 



Ell |Vo 
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where a — 1, . . . , N indexes the state variables of numerical solution vector U and the reference solution V, 
and i ~ 1, . . . , / indexes the grid cells of the entire computational domain. For test problems with smooth 
solutions, we will also use the relative maximum error defined by 

P V" maxj=i,...,j|Ua— Vail 

""Loo-I. max..,...,|V„,| ' ^'^^ 

Quite often the reference solution is defined on a grid with higher resolution than that of the numerical 
solution. In that case, we first coarsen the reference solution to the resolution of the numerical solution. 



4.2. Radiation Tests 

4-.2.1. Su-Olson Test 



Su fc Olson I (|l996f ) developed a one-dimensional Marshak wave test, to check the accuracy of the scheme 
and the correctness of the implementation of the time-dependent non-equilibrium gray radiation diffusion 
model. In this test, radiation propagates through a cold medium that is initially in absence of radiation. The 
equations are linearized by the choice of the specific heat of the material Cy = AaT^ as well as by setting 
the Rosseland and Planck opacities to the same uniform and time-independent constant kji = Kp = k. The 
cold medium is defined on a half-space of the slab geometry < x < oo. At the boundary on the left, 
a radiative source is specified, creating an incident radiation flux of i^"' — aT^^^ where Tin = 1 keV. As 
time progresses, the radiation diffuses through the initially cold medium and b y energy exchange between 



radiation and matter, the material temperature will rise. In ISu fc Olson I (|l996r ). a semi-analytical solution 
is derived for the time evolution of the radiation energy and material temperature. We will use this solution 
for this verification test. 

For convenience, we locate the right boundary at the finite distance a; = 5 cm and impose a zero 
incoming radiation flux on that boundary. We decompose the computational domain in 6 grid blocks at 
the base level with 10 cells per block. Between a; = 5/6 cm and a; = 5/3 cm, the domain is refined by one 
level of AMR. During the time evolution, the radiation diffuses to the right through the resolution changes. 
The system is time evolved with the implicit radiation diffusion solver by using a preconditioned conjugate 
gradient method till the final time 0.02 ns. The solver steps through a series of fixed time steps of 5 x 10~^ ns 
and we use a Crank-Nicolson approach to achieve second order accurate time-integration. Note that this is 
possible because coefficients of the matrix to be solved are not time dependent. The computed radiation 
and material temperatures at the final time are shown in Figure [1] and agree well with the semi-analytical 
solution. 

Figure shows the relative LI error of the radiation and material teperatures versus increasing grid 
resolution of the base level grid. We did not use the semi-analytical solution as the reference, since it i s 



difficult to get an accurate enough solution with the quadrature method as mentioned bv lSu fc Olson I (jl996l ). 
Instead, we use a very high resolution (1920 cells) numerical reference solution obtained with the CRASH 
code. Four different base level resolutions with 60, 120, 240, and 480 cells are used to demonstrate the second 
order convergence. The time step is proportional to the cell size Ax. 



- 20 - 



4-2.2. Lowrie's N on- equilibrium Radhydro Solutions 



Lowrie fc Edwards I (|2008l ) designed several shock tube problems for the non-equilibrium gray radiation 
diffusion coupled to the hydrodynamic equations that can be used for code verification. These solutions are 
planar radiative shock waves where the material and radiation temperatures are out of equilibrium near the 
shock, but far from the shock the flow is assumed to be in radiative equilibrium. Depending on the Mach 
number of the pre-shock state, a wide range of shock behavior can occur. For the CRASH test suite, we 



selected a few of the semi- analytic solutions from lLowrie fc Edwards I (|2008l ). In this section we will describe 



the Mach 1.05 flow with uniform opacities as an example. Here the shock is smoothed out by the energy 
exchange with the diffusive radiation. Another more challenging Mach 5 problem with non-uniform opacities 
will be described in Section 14.3.31 

The Mach 1.05 test is performed on a 2D non- uniform grid. The initial condition is taken to be the 
same as the original steady state reference solution. Since the system of equations is Galilean invariant, we 
can add an additional velocity -1.05 so that the velocity on the left boundary is zero while the smoothed 
shock will now move to the left. This new initial condition as well as the velocity vector are rotated by 
tan^^(l/2) w 26.56°. This means that there is a translational symmetry in the (—1,2) direction of the 
xy-plane as shown in Figure [31 The computational domain is —0.12 < x < 0.12 by —0.02 < y < 0.02 
decomposed in 3 x 3 grid blocks of 24 x 4 cells each. We apply one level of refinement inside the region 
-0.04 < X < 0.04 by -0.02/3 < y < 0.02/3. The initial smoothed shock starts at the right boundary of 
the refined grid and we time evolve the solution till it reaches the resolution change on the left as shown 
in Figure [3l For the boundary conditions in the x direction we use zero radiation infiux conditions for the 
radiation field, while zero gradient is applied to the remaining state variables. On the y boundaries, we apply 
a sheared zero gradient in the (—1,2) direction for all variables. 

The hydrodynamic equations are time evolved with the HLLE scheme with a CFL number 0.8. We use 
the generalized Koren limiter with /? = 3/2 for the slope reconstruction. For the implicit radiation diffusion 
solver, we use the GMRES iterative solver in combination with a BILU preconditioner. The specific heat is 
time dependent since it depends on the density, therefore the implicit scheme is only first order accurate in 
time. To enable second order grid convergence for this smooth test problem, we compensate this by reducing 
CFL number proportional to the grid cell size, in other words At cx Ax^, so that second order accuracy with 
respect to Aa; can be achieved. We increase the spatial resolution by each time doubling the number of grid 
blocks at the base level in both the x and y directions. 

The convergence of the numerically obtained material and radiation temperatures along the y = cut at 
the final time t = 0.07 is shown in Figure ID The solid, dotted, and dashed lines correspond to the solutions 
with the 3 X 3, 6 X 6, and 12 x 12 base level grid blocks, respectively. The advected semi-analytical reference 
solution is shown as a blue line for comparison. 

To assess the order of accuracy, the grid convergence is shown in Figure [5] for the three resolutions. The 
relative LI error is calculated using the density, velocity components, and both the material and radiation 
temperatures. We obtain second order convergence for both the conservative as well as the non-conservative 
(using the pressure equation instead of the total energy) hydrodynamic schemes. The latter scheme can be 
used because in the Mach 1.05 test the hydro shock is smoothed out by the interaction with the radiation. 
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4-2.3. Double Light Front 

As a test for the multigroup radiation diffusion model we developed a double light front test problem. 
This test is used to verify the implementation of both the group diffusion and flux limiters. At the light 
front the discontinuity in the radiation field switches on the flux limiter. This limiter is used to correct the 
radiation propagation speed in the optically thin free streaming regime. With the light front test we can 
then check that we obtain the speed of light propagation of the front and that the front maintains as much 
as possible the initial discontinuity. 

This test is constructed as follows: We use a ID computational domain of the size of 1 m in the a;- 
direction. On this domain, we initialize the two radiation group energy densities Eg (g = 1,2) with a very 
small, positive number to avoid division by zero in the flux limiter. Also the Rosseland mean opacities are 
set to a small number corresponding to strong radiation diffusion, while the Planck mean opacities are set to 
zero corresponding to an optically thin medium. The radiation energy density of the first group enters from 
the left boundary by applying a fixed boundary condition with value one in arbitrary units. On the right 
boundary this group is extrapolated with zero gradient. Note that these are the proper boundary conditions 
in the free-streaming limit and not the diffusive flux boundary conditions described in Section 13.41 The 
second radiation group enters from the right boundary with density one, and it is extrapolated with zero 
gradient at the left boundary. We time evolve both groups for 0.5 m/c seconds. The analytic solution will 
then be two discontinuities that have reached x = 0.5 m, since both fronts propagate with the speed of light 
c. 

The computational domain is non-uniform. In the coarsest resolution there are 10 grid blocks of 4 cells 
each at the base level. Inside the regions 0.1 < a; < 0.2 and 0.8 < a; < 0.9, we use one level of refinement. 
The total time evolution is divided into 400 fixed time steps. We use GMRES for the radiation diffusion 
solver in combination with a BILU preconditioner. For the grid convergence we reduce the fixed time step 
quadratically with the grid resolution. This time step reduction mimics second order discretization in time. 
In Figureini the two group energy densities are shown for the base level grid resolutions 40, 80, 160, and 320. 
Clearly, with increasing number of cells, the solution converges towards the reference discontinuous fronts at 
X = 0.5. 

In Figure [71 the grid convergence is shown for the four resolutions. The relative LI error is calculated 
using both radiation group energy densities and compared to the analytical reference solution with the 



discontinuities at x = 0.5. In lGittings et al. I (|2008l ). it was stated that for a second order difference scheme 
the convergence rate for a contact discontinuity is 2/3. Indeed, we find this type of convergence rate, due to 
the numerical diffusion of the discontinuities, for the light front test. We have also performed the tests in 
the y and z directions to further verify the implementation. 



4.2.4. Relaxation of Radiation Energy Test 

This test is designed to check the relaxation rate between the material and the radiation. The energy 
exchange between the material and radiation groups can be written as 

dT ^ 

= Y.^9{E,-Bg), (90) 
9=1 

BE 

= Og{Bg-Eg). (91) 
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For a single radiation group, an analytic expression can be found to describe the relaxation in time. However, 
for arbitrary number of groups, a time dependent analytic solution is less obvious, except for some rather 
artificial cases. Here we make the assumption of an extremely large value of the specific heat Cv to make 
the analysis more tractable. In that case, the material temperature is time independent, so that Bg is 
likewise time independent. The solution is then Eg = Bg{\ — e^'^'*) assuming Eg{t = 0) = initially. At 
time t = l/<Jg, the group radiation energy density is Eg = Bg(l — 1/e). Note that this test only needs one 
computational mesh cell in the spatial domain. We set T = IkeV and the resulting Planckian spectrum, 
defined by Bg , is depicted by the dotted line in the left panel of Figure [H We use 80 groups logarithmically 
distributed over the photon energy domain in the range of 0.1 eV to 20keV. The computed Eg at time 
t = \/(Jg are shown as + points. For the simulation we used the GMRES iterative solver and the Crank- 
Nicolson scheme. To assess the error, we repeated the test with time steps of 1/20, 1/40, and 1/80 of the 
simulation time. The second order convergence rate is demonstrated in the right panel of Figure [S] 



4.3. Heat Conduction Tests 

4-. 3.1. Uniform Heat Conduction in rz- geometry 

This test is designed to verify the implicit heat conduction solver in rz-geometry. It tests the time 
evolution of the temperature profile using uniform and time independent heat conductivity. In rz-geometry, 
the equation of the electron temperature for purely heat conductive plasma follows 

Cy,^ = i A (rC,^) + ^ (c,^) . (92) 
dt r dr \ dr ) dz \ dz ) 

We will set the electron specific heat Cve = 1 and assume the electron conductivity Ce to be constant. In 
that case, a solution can be written as a product of a Gaussian profile in the z-direction and an elevated 
Bessel function Jo in the r-direction ( Arfken 19851 ): 



Te = T„,i„ + To-=i==e-3!^ Jo(He"''''^'=*, (93) 

where h w 3.8317 is the first root of the derivative of Jo{r). We select the following values for the input 
parameters: Ti^in = 3, Tq = 10, and = 0.1 in dimensionless units. 

The computational domain is — 3 < z < 3 and < r < 1 discretized with 3x3 grid blocks of 30 x 30 
cells each. In the region — 1 < z < 1 and 1/3 < r < 2/3, we apply one level of mesh refinement. We 
impose a symmetry condition for the electron temperature on the axis. On all other boundaries the electron 
temperature is fixed to the time dependent reference solution. We time evolve this heat conduction problem 
with a preconditioned conjugate gradient method from time t = 1 to the final time at t = 1.5. The Crank- 
Nicolson approach is used to achieve second order accurate time integration. 

The initial and final solutions for the electron temperature are shown in Figure [5] in color contour in 
the 7'z-plane. The heat conduction has diffused the temperature in time to a more uniform state. The 
black line indicates the region in which the mesh refinement was applied. The relative maximum error of 
the numerically obtained electron temperature versus the analytical solution is shown in Figure 1101 Here 
we used the non- uniform grid with base resolutions of 90^, 180^, 360^, and 720^ cells and set the time step 
proportional to the cell size to demonstrate a second order convergence rate. 
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.3.2. Reinicke Meyer-ter Vehn Test 



The iReinicke fc Mever-ter-Vehn I ()1991l) problem tests both the hydrodynamics and the heat conduc- 
tion implementation. This test generalizes the well-known Sedov- Taylor strong point explosion in single 
temperature hydrodynamics by including the heat conduction. The heat conductivity is parameterized as a 
non- linear function of the density and material temperature: Ce = p'^T^. We select the spherically symmet- 



ric self-similar solution of IReinicke fc Mever-ter-Vehn I (|1991I ) with coefficients a — —2 and b — 13/2 and the 



adiabatic index is 7 — 5/4. This solution produces, similar to the Sedov-Taylor blast- wave, an expanding 
shock front through an ambient medium. However, at very high temperatures, thermal heat conduction 
dominates the fluid flows, so that a thermal front precedes the shock front. With the selected parameters, 
the heat front is always at twice the distance from the origin of the explosion as is the shock front. 

We perform the test in rz-geometry. The computational domain is divided in 200 x 200 cells. The 
boundary conditions along the r and z axes are reflective. The two other boundaries, away from explosion, 
are prescribed by the self-similar solution. The time evolution is numerically performed as follows: For the 
hydrodynamics we use the HLLE scheme and the generalized Koren limiter with /3 = 2 as the slope limiter. 
The CFL number is set to 0.8. The heat conduction is solved implicitly with the preconditioned conjugate 
gradient method. The test is initialized with the spherical self-similar solution with the shock front located 
at the spherical radius 0.225 and the heat front is at 0.45. The simulation is stopped once the shock front 
has reached 0.45 and the heat front is at 0.9. 

A ID slice along the r-axis of the solution a t the final time is shown in Figure [TT] We normalize 



the output similar to IReinicke fc Mever-ter-Vehn I (|1991I) : The temperature is normalized by the central 



temperature, while the density and radial velocity are normalized by their values of the post-shock state at 
the shock. The numerical solution obtained by the CRASH code is shown as -I- symbols and is close to the 
self-similar reference solution, shown as solid lines. Note that the temperature is smooth due to the heat 
conduction, except for the discontinuous derivative at the heat front. The wiggle at r = 0.3 in the density 
and radial velocity is a due to the diffusion of the analytical shock discontinuity in the initial condition 
during the first few time steps. In the left panel of Figure 1121 the spherical expansion of temperature at 
the final time is shown. Clearly, the Cartesian grid with the rz-geometry does not significantly distort the 
spherical symmetry of the solution. The spatial distribution of the error in the temperature is shown in the 
right panel. The errors are largest at the discontinuities of the shock and heat fronts as expected. 

A grid convergence study is performed with resolutions of 200^, 400^, and 800^ cells. The relative LI 
error in Figure |T3| is calculated using the density, velocity components, and the material temperature. The 
convergence rate is first order due to the shock and heat front. 



4.3.3. Heat Conduction Version of Lowrie's Test 



Any of the verification tests for non-equilibrium gray-diffusion coupled to the single temperature hy- 
drodynamics can be reworked as a test for the hydrodynamic equations for the ions coupled to the electron 
pressure equation with electron heat conduction and energy exchange betwee n the electrons and ions. As an 
example, we will transform one of the non-equilibrium gray diffusion tests of iLowrie fc Edwards I (|2008l ) to 
verify the heat conduction implementation. 
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The electron energy density equation (j20l) without the radiation interaction can be written as 

^ + V • [E,u] + peV • U = V • [CeVTe] + a,e(T, - Te), (94) 

where the heat conduction and energy exchange terms on the right hand side depend on the gradients and 
differences of the temperatures. The equation for the gray radiation energy density (|16p on the other hand 
depends on the gradients and differences of energy densities. By defining the radiation temperature Tr by 
Er = aT^ and using the definition of the Planckian B = aT^, we can rewrite the energy density equation 
for the radiation as 

OE 1 

—^+V-\ErU] + -ErV-U = V ■[DrVEr\+CKp{aT'^ - Er) 

at 3 

= V • [A-VT^] + cKp(T - Tr), (95) 

where Dr = Dr-iaTr and CKp = cKpa{T^ + Tr)(T + Tr) arc the new coefficients that appear due to this 
transformation. The equations (j94p and (j95p are now of the same form. To translate a gray diffusion test to 
a heat conduction test, we reinterpret Dr as the heat conductivity Cg and CKp as the relaxation coefficient 
aie in the ion-electron energy exchange. In addition, the material temperature T and radiation temperature 
Tr have to be reinterpreted as the ion temperature Ti and electron temperature T^, respectively. Note that 
we also have to relate the electron pressure and internal energy by — Ee/3 similar to the radiation field 
corresponding to 7e = 4/3, and let the electron internal energy and electron specific heat depend on the 
electron temperature as Eg = aT^ and Cve — ^clT^i respectively. 



As an example, we transform the Mach 5 non-equilibrium gray diffusion shock tube problem of lLowrie fc Edwards 



( 20081 ). It uses non-uniform opacities that depend on the density and temperature defined by D 



0.0175(7r)^/Vp and CKp = lO^/Z^r- The above described procedure is used to translate this problem 
to an electron heat conduction test with energy exchange between the electron and ions. The heat conduc- 
tivity for this test is Ce = ■iaT^{).Qnb{'-^Ti)'^ /"^ / p and the relaxation coefficient between the electron and ions 
is a,e = a{Tf + T2)(r, + T,)AaT^lQ'' / C^. 

We perform this Mach 5 heat conduction test o n a 2D non-uniform grid. F or the initial condition, the 



ID semi-analytical steady state reference solution of lLowrie fc Edwards I (j2008l ) is used. There is a Mach 5 



pre-shock flow on the left side of the tube resulting in an embedded hydro shock as well as a steep thermal 
front (a look at Figure [14] will help to understand this shock tube problem.) We add an additional velocity 
of Mach —5 so that the pre-shock velocity is zero and the shock is no longer steady, but instead will move 
to the left with a velocity —5 (in units in which the pre-shock speed of sound is 1). The problem is rotated 
anti-clockwise on the grid by tan~^(l/2). The translational symmetry is now in the (—1,2) direction in 
the xy-plane similar to the Mach 1.05 shock tube problem described in Section [4.2.21 The computational 
domain is -0.0384 < x < 0.0384 by -0.0048 < y < 0.0048. Inside the area -0.0128 < x < 0.0128 and 
—0.0016 < y < 0.0016, we apply one level of refinement. This refinement is set up such that both the heat 
front as well as the shock front will propagate through the resolution change on the left (from fine to coarse) 
and right (from coarse to fine), respectively. For the boundary conditions in the x-direction, we fix the state 
on the right side with the semi-analytical solution, but for the left side we use zero gradient. On the y 
boundaries, we apply a sheared zero gradient in the (—1,2) direction. 

For the evolution till the final time t — 0.0025, we use the HLLE scheme together with the generalized 
Koren limiter with /3 = 3/2 to solve the hydrodynamic equations. The CFL number is set to 0.8. The heat 
conduction and energy exchange between electrons and ions are solved implicitly with the backward Euler 
scheme using the GMRES iterative solver in combination with a BILU prcconditioner. 
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In Figure 1141 the electron (right panel) and ion (left panel) temperatures are shown at the final time 
along the x-axis. The semi-analytical reference solution is shown as a blue line, while the numerical solution 
is shown with + symbols for a simulation with 192 x 24 cells at the base level in the x and y direction. The 
hydro shock is located near x w 0.0085 and shows up in the ion temperature as a jump in the temperature, 
followed directly behind the shock by a strong relaxation due to the energy exchange between the ions and 
electrons. The electron temperature stays smooth due to strong heat conduction. The heat front is seen 
with a steep foot at x ~ —0.022. This fro n t corr esponds to the radiative precursor in the non-equilibirum 



gray diffusion tests of iLowrie fc Edwards I (|2008l ). We repeated the test with four different resolutions at 
the base level: 192 x 24, 384 x 48, 768 x 96, and 1536 x 192 cells in the x and y direction. The insets in 
both panels of Figure [Ml show the four resolutions as solid, dotted, dashed, dashed-dotted fines, respectively. 
In the left panel, the zoom-in shows the convergence of the ion temperature towards the embedded hydro 
shock and the temperature relaxation. In the right panel, the blow-up shows the convergence towards the 
reference precursor front. Note that no spurious oscillations appear near the shock or near the precursor. 

Due to the discontinuity in both the shock and heat precursor, the convergence rate can be at most first 
order. Indeed, in Figure [15] the relative LI error shows first order accuracy. The error is calculated using 
all the density, velocity components and both temperatures. Note that the spike in the ion temperature is 
spatially so small that a huge number of grid cells are needed to get a fully resolved shock and relaxation 
state. 



4.4. Full System Tests 

The CRASH test repository contains a range of full system configurations to be used for validation with 
future laboratory experiments. In Figure [TBI we show the configuration of a 3D elliptic nozzle through which 
a fast shock of the order of 150km/s will be launched, which is still significantly slower than the speed of 
light. The shock wave is produced by a 1.1ns laser pulse from the left with 4kJ of energy irradiating a 
20 /im thick Beryllium disk, initially located at x — 0. A layer of gold is glued to the plastic tube to protect 
the outside of the tube from the laser-driven shock. The plastic (polyimidc) tube is circular for x < 500 /im 
with a radius of 600 /im. Beyond x ~ 750 /im the tube is elliptic by fiattening the tube in the y-direction by 
a factor 2. 

The fi rst part of the simula tion is performed with the 2D, Lagrangian, radiation hydrodynamic code 



HYADES (jLarsen fc Lane Ill994[ ) to time advance the laser energy deposition and the response of the system 
until the end of the laser pulse at 1.1 ns. This laser pulse first shocks and then accelerates the Beryllium to 
the right. After 1.1 ns the output of HYADES is used as an initial condition of the CRASH code. 

This simulation is performed for a two temperature, electron and ion, plasma. For the radiation, we 
use the flux limited diffusion approximation with 30 groups. The photon energy is in the range of 0.1 eV to 
20keV, logarithmically distributed over the groups. Due to the symmetry in the problem we only simulate 
one quadrant (jj > and z > 0), with reflective boundary conditions at y = and z — 0. At all other 
boundaries we use an extrapolation with zero gradient for the plasma and a zero incoming flux boundary for 
the radiation. The domain size is [—150,3900] x [0,900] x [0,900] microns for the x, y, z coordinates. The 
base level grid consists of 120 x 20 x 20 blocks of 4 x 4 x 4 mesh cells. One level of dynamic mesh refinement 
is used at material interfaces and the shock front. Overall, the effective resolution is 960 x 160 x 160 cells 
and there are approximately 4.5 million finite volume cells. The hydrodynamic equations are solved with 
the HLLE scheme with a CFL number 0.8 together with the generalized Koren limiter with /? — 3/2. The 
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diffusion and energy exchange of the radiation groups as well as the heat conduction are solved with the 
decoupled implicit scheme using a Bi-CGSTAB iterative solver. The simulation from 1.1ns to 13 ns physical 
time took 1 hour and 55 minutes on 480 cores of the FLUX supercomputer at the University of Michigan. 

In Figure [TTl we show the shock structure at 13 ns. The accelerated Beryllium compresses the Xenon 
directly to the right of the interface, which is seen as a high density plasma near x — 1700/im in the top 
right panel of Figure [TTl This drives a primary shock and the velocity jump at a; ~ 1700 /^m is seen in 
the middle left panel. Behind the shock front, the ions are heated as depicted from the middle right panel, 
followed directly behind the shock by a cooling due to the energy exchange between the ions and electrons. 
Early on, the electron heating produces ionization and the emission of radiation, and the radiation in turn 
heat and ionizes the material ahead of the primary shock. The radiation temperature, measuring the total 
radiation energy density, is shown in the bottom left panel. The photons will interact again with the matter, 
someti mes after traveliiig sorn e distance. This is the source of the wall shock seen ahead of the primary 
shock ( Doss et al. 20091 2011 ): photons traveling ahead of the shock interact with the plastic wall, heat it, 
and this in turn drives a shock off the wall into the Xenon. The ablation of the plastic is depicted in the top 
left panel as a radially inward moving polyimide (in green color) near and even ahead of the primary shock. 
The compressed Xenon due to the plastic ablation is seen in the top right panel as a faint density feature 
that is ahead of the primary shock front, between x — 1700 /im and x — 2000 /im. The interaction between 
the photons and matter is also seen by the radiative precursor to the right of the radiative shock elevating 
the electron temperature ahead of the shock in the bottom right p anel. This is d ue to the strong coupling 
between the electrons and radiation field. The reader is referred to iDrake et al. I (j201ll ) for more details on 
radiative effects in radiative shock tubes. 



4.5. Parallel Performance 



We present parallel scaling studies on the Pleiades supercomputer at NASA Ames. This computer is 
an SGI ICE cluster connected with infiniband. Figure [18] shows the strong scaling for a problem size that 
is independent of the number of processors. This 3D simulation is a circular tube version of the full system 
test described in Section 14.41 It uses five materials, 30 radiation groups, and separate electron and ion 
temperatures. The grid contains 80 x 8 x 8 blocks of 4 x 4 x 4 cells each at the base level and in addition 
two time dependent refinement levels. There are overall approximately 2.6 million cells in this problem. 
We use lookup tables for the EOS and opacities, so that the computational time for that is negligible. For 
the hydrodynamic equations, we use the HLLE scheme together with the generalized Koren limiter with 
P = 3/2. The radiation diffusion, electron heat conduction and energy exchange terms are solved implicitly 
with the decoupled scheme, uing the Bi-CGSTAB iterative solver. This simulation is performed for 20 time 
steps for the number of cores varying from 128 to 2048, but excludes file I/O to measure the performance of 
the implicit solver. Up to 1024 cores, we get good scaling. However for more cores we observe saturation in 
the performance. 



SUMMARY 



We have extended the BATS-R-US code |Powell et al. I Il999t iToth et al. I I2OI0I ) with a new radia- 
tion transfer and heat conduction library. This new combination together with the equation-of-state and 
multigroup opacity solver is called the CRASH code. This code uses the recently developed parallel Block 



Adaptive Tree Library (BATL, see iToth et al. I (|2010l )) to enable highly resolved radiation hydrodynamic 
solutions. The implemented radiation hydrodynamic schemes solve for the gray or niultigroup radiation 
diffusion models in the flux limited diffusion approximation. 



In high energy density plasmas, the electrons are most of the time strongly coupled to the ions by 
collisions. An important exception is at hydrodynamic shocks, where the ions are heated by the shock 
wave and the electrons and ions are out of temperature equilibrium. Since radiative shocks are the main 
application for CRASH, we have implemented a separate electron pressure equation with the electron thermal 
heat conduction. For the electron heat conduction, we have added the option of a flux limiter to limit the 
thermal flux with the free-streaming heat flux. 

The multi-material radiation hydrodynamic equations are solved with an operator split method that 
consists of three substeps: (1) solving the hydrodynamic equations with standard finite volume shock- 
capturing schemes, (2) the linear advection of the radiation in frequency-logarithm space, and (3) the implicit 
solution of the radiation, heat conduction, and energy exchanges. For the implicit solver, standard Krylov 
solvers are used together with a Block Incomplete Lower-Upper decomposition (BILU) preconditioner. This 
preconditioner scales well up to 500 or 1000 processors. For future work, we may explore for the implicit 
multigroup diffusion a multi-level preconditioner to better scale the radiation solver beyond 1000 processors. 

We have presented a suite of verification tests that benchmark the performance. These tests verify the 
correctness and accuracy of the implementation for the gray and multigroup radiation diffusion algorithm 
and the heat conduction in ID, 2D, and 3D slab and 2D rz geometry. To demonstrate the full capability of 
the implementation, we have presented a 3D multi-material simulation of a radiative shock wave propagating 
through an elliptic nozzle. This configuration will be used in future validation studies. 

Since this radiation transfer library is an extension of the BATS-R-US code, the implementation is 
readily available for MHD simulations as well. This allows for validation studies of the radiation MHD 
implementation using laboratory-astrophysics experiments or for the simulations of astrophysical plasmas. 

This work was funded by the Predictive Sciences Academic Alliances Program in DOE/NNSA-ASC 
via grant DEFC52-08NA28616 and by the University of Michigan. The simulations were performed on the 
NASA Advanced Supercomputing system Pleiades. 



A. DISCRETIZATION OF THE DIFFUSION OPERATOR AT RESOLUTION 

CHANGES 

In Sections 13.3.11 and 13.3.21 the diffusion operator is discretized on a uniform mesh with a standard 
finite volume method in combination with a central difference approximation for the gradient in the flux 
calculation as in equation (|68|) . The diffusion coefficient that is needed at the face is obtained by simple 
arithmetic averaging of the left and right cell center diffusion coefficients. The generalization to resolution 
changes as in Figure [19] is less straightforward. In the following, we will denote the fine cell centers by a and 
b, the coarse cell center by c. The flux densities at the resolution changes in the direction orthogonal to the 
interface are denoted by Fi and F2 at the flue faces, and F3 at the coarse face. 



In I Edwards I (Il996[ ). a strategy was developed to discretize the diffusion operator on an adaptive mesh 
in the context of reservoir simulations. The main ingredients of the method are (1) require the continuity 
of the flux at the resolution change in the strong sense, i.e. Fi = F2 = ^3, and (2) discretize the gradient 
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in the diffusion flux by a one-sided difference. An expression was found for the diffusion flux F = —DVE 
in which the diffusion coeffi c ient is replaced by a weighted harmonic average of the cell centered values Da, 
Db, Dc- In [cittings et al. I ( 20081) , it was argued that this discretization does not properly propagate the 
self-similar Marshak waves of the radiation diffusion model, unless the cell centered diffusion coefficients are 
calculated on a common face temperature. 

In the code discussed in this paper, we follow a different approach that replaces the harmonic average of 
the diffusion coefficient in lEdwards I |l996l) by an arithmetic average and obtain for the flux densities normal 
to the resolution change interface 



Fi — F2 — 



2D 

"3A^ 



[E, - {Fa + Eb)/2] , 



where Aa; is the fine cell size and the diffusion coefficient at the face is averaged as 

D3 = [D, + {Da + Db)] /3. 



(Al) 



(A2) 



We demonstrated with verification tests including those discussed above that this change produces properly 
propagating radiative precursor and shock fronts. Generalizations to ID and 3D are straightforward. 



B. RZ-GEOMETRY 

Incorporating the rz-geometry in a finite volume formulation is as follows: the radial cell face area and 
the cell volume must be made proportional to the distance r from the symmetry axis. In addition, the r 
component of the momentum equation (|18p is modified as 

—g^+V-[puur + r{p + pr)] = — - — , (Bl) 

where f is the unit vector in the r direction and Ur = u-r. This correction refiects the fact that the pressure 
term is a gradient, not a divergence. 
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Table 1. Quantities stored in the EOS tables as a function of logTe [eV] and logria [m ^] . 



quantity 


stored quantity 


units 


total pressure p 


p/na 


eV 


total internal energy density E 


E/na 


eV 


electron pressure Pe 


Pe/ria 


eV 


electron internal energy density E,, 


Ee/Ua 


eV 


specific heat Cy 


Cv/{nakB) 




electron specific heat Cve 


Cve/inaks) 




speed of sound gamma 


7s 




electron speed of sound gamma 






inverse of ion-electron interaction time 




s-i 


electron conductivity 


C'e 


Jm-^s-iR-i 


mean ionization 


Z 




moan s(|uar(> ionization 
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Table 2. Quantities stored in the opacity tables as a function of logp [kg m "^j and logTg [eV]. 



quantity 


symbol 


units 


specific group Rosseland mean opacities 




kg^^ m^ 


specific group Planck mean opacities 


Kpg/p 


kg~^ m^ 
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Fig. 1. — The material {T-^g^^) and radiation (T-^g^f^) temperature solution of the ISu fc Olson I (Il996l) non- 
equilibrium Marshak radiation diffusion problem obtai ned with the CRAS H code on a non-uniform grid. 
The reference temperatures of the analytical method of ISu fc Olson (1996) are shown as lines. 
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Fig. 2. — The relative LI error for the Su-Olson test on a non-uniform grid. 
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Fig. 3. — Rotated shock tube te st on a 2D AMR grid based on the Mach number 1.05 non-equilibrium gray 
radiation hydrodynamic test in iLowrie fc Edwards I (|2008l ). Shown is the radiation temperature in color 
contour at the initial (top panel) and final (bottom panel) times. The black crosses indicate the cell centers. 
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Fig. 4. — The material (left panel) and radiation (right panel) temperatures for the Mach 1.05 radiative 
shock tube problem at the final time are shown in the x-direction. The solid, dotted, and dashed lines 
correspond to three different grid re solutions, respectively. The blue line is the semi-analytical reference 



solution of Lowrie fc Edwards ( 20081) . 
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Fig. 5. — The relative LI error for the Mach 1.05 non-equihbrium radiation diffusion test on a non-uniform 
grid. Both the non-conservative as well as the conservative hydrodynamic schemes are tested. 
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Fig. 6. Solutions for the ID double light front test for 4 different non-uniform grid resolutions. The 
radiation energy for group 1 (left panel) enters from the left boundary, for group 2 (right panel) it enters from 
the right boundary. The symbols for base resolution 80 shows one level of grid refinement for 0.1 < x < 0.2 
and 0.8 < a; < 0.9. 




Fig. 7. — Relative LI error for the double light front test on a non- uniform grid. The test is performed for 
the X, y, and z directions. 
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Fig. 8. — The relaxation of radiation energy test for 80 groups. Left panel is for the time independent 
spectrum Bg (dotted line) and the group radiation energy solution Eg at time l/cTg (+ points) versus the 
photon energies after 80 time steps. The analytical reference solution is shown as a solid line. Right panel 
shows the relative maximum error for 20, 40, and 80 time steps demonstrating second order convergence 
rate. 
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Fig. 9. — The electron temperature for the uniform heat conduction test on a non-uniform grid in rz- 
geometry. The top panel shows the electron temperature in the initial condition while the bottom panel 
is the electron temperature at the final time. The black box indicates the region within which the grid is 
refined by one level. 
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Fig. 10. The relative maximum error for the uniform heat conduction test on a non-uniform grid in the 
r^;-geometry. 




Fig. 11. — Density (top panel), temperature (middle panel), and radial velocity (bottom panel) along the 

z = cut for the Reinicke Meyer-ter Vehn test in rz-geometry. The numerical solution (+ symbols) is at 
the final time compared to the self-similar analytical reference solution (solid lines). 
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Fig. 12. — The temperature (left panel) and temperature error compared to the reference solution (right 
panel) for the Reinicke Meyer-ter Vehn test in rz-geometry. 
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Fig. 14. — Mach 5 shock tube problem of iLowrie &: Edwards I (|2008t l transformed to a non-uniform heat 
conduction and ion-electron coUison frequency test and rotated on a 2D non-uniform grid. The ion (left 
panel) and electron (right panel) temperatures at the final time are shown in the x-direction. The blue line 
is the reference solution. In the left panel, the grid convergence near the shock is shown in the inset. In the 
right panel, a blow-up of the grid convergence to the reference heat front is shown. 
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Fig. 15. — Relative LI error for the Mach 5 non-equilibrium heat conduction test on a non-uniform grid. 
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Fig. 16. — The geometry of the 3D eUiptic nozzle experiment after 1.1 ns, consisting of 5 materials: Beryllium 
(blue), Xenon (black), polyimide (green), gold (yellow), and acrylic (red) in both panels. The radius of the 
inside of the polyimide tube is 600 /zm in the y = plane (left panel). In the z = plane (right panel), the 
radius of the inner tube is 600 /xm for x < 500 /im, but shrinks to 300 /im beyond x ~ 750 /im. The lines 
represent the mesh refinement at material interfaces and shock fronts. 
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Fig. 17. — Simulated radiative shock structure at 13ns in a 3D elliptic nozzle consisting of the 5 materials 
indicated in Figure 1161 The plots show in the xy-plane in color contour the variables indicated in the plot 
title. The primary shock is at x « 1700. 



- 50 - 




Fig. 18. — Strong scaling of the CRASH code, running a 3D CRASH application with 5 material level sets, 
electron and ion temperature, 30 radiation groups, and two levels of time dependent mesh refinement. 
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Fig. 19. — Cell and face centers at the adaptive interface in 2D. 



